Accurate Parallel Integration of Large Sparse Systems of Di erential Equations

We describe a MIMD parallel code to solve a general class of ordinary di erential equations, with particular emphasis on the large, sparse systems arising from space discretization of systems of parabolic partial di erential equations. The main goals of this work are sharp bounds on the accuracy of the computed solution and exibility of the software. We discuss the sources of error in solving di erential equations, and the resulting constraints on time steps. We also discuss the theory of a posteriori error analysis for the Galerkin nite element methods, and its implementation in error control and estimation. The software is designed in a matrix-free fashion, which enables the solver to e ectively tackle large sparse systems with minimal memory consumption and an easy and natural transition to MIMD (distributed memory) parallelism. In addition, there is no need for the choice of a particular representation of a sparse matrix. All memory is dynamically allocated, with a new expandable array object used for archiving results. The implicit solution of the discrete equations is carried out by replaceable modules: the nonlinear solver module may be a full Newton scheme or a quasi-Newton; these in turn are implemented with a linear solver, for which we have used both a direct solver and QMR, an iterative (Krylov space) method. Three example computations are presented: the Lorenz system, which has dimension three and the discretized versions of the (partial-di erential) bistable equation in one and two dimensions. The Lorenz system demonstrates the quality of the error estimation. The discretized bistable examples provide large sparse systems, and our precise error estimation shows, contrary to standard error estimates, that reliable computation is possible for large times.